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»' I We introduce an urn model which describes spatial separation of sand. In this dynamical model, 

in a certain range of parameters spontaneous symmetry breaking takes place and equipartitioning 
^ ' of sand into two compartments is broken. The steady-state equation for an order parameter, a 

O . critical line, and the tricritical point on the phase diagram are found exactly. Master equation and 

■ the first-passage problem for the model are solved numerically and the results are used to locate 

first-order transitions. Exponential divergence of a certain characteristic time shows that the model 
can also exhibit very strong metastability. In certain cases characteristic time diverges as A*'^ , where 
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N is the number of balls and z — (critical line), | (tricritical point), or i (limits of stability). 



o 
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I. INTRODUCTION 

c/5 I Recently, granular systems are intensively studied using methods of statistical physics. The main motivation of this 
research is to get basic understanding of these technologically important materials whose behaviour, however, very 
often appears elusive and mysterious ||^,|^. 
Q , One of the classical experiments in this field concerns the spatial separation of shaken sand |^ . In this experiment 
' . ' one uses a box which is divided into two equal compartments by a wall which has a narrow horizontal slit at a certain 
^ height. When the box is filled with a sand and shaken one observes that under certain conditions (e.g., frequency and 
Q ■ amplitude of shaking) one compartment of a box is preferentially filled with sand. Theoretical arguments to understand 
' how the symmetry breaking arises in this experiment, were given by Eggers In his approach, Eggers studied a 
continuous model based on partial differential equations which takes into account certain particular characteristics 
of shaken granular systems Namely the fact that effective temperature of a granular system decreases when 
^ the density of particles increases. This at first counterintuitive feature of granular systems is now understood as a 
^\ result of inelastic collisions between particles which effectively cool the system This feature is also responsible 
0^ for the symmetry breaking in the experiment with shaken sand. Indeed, with nonequal distribution of sand the 
\l temperature in the compartment with majority of particles is substantially reduced. Consequently, particles from this 
j ' compartment are less likely to cross the slit. Eggers has shown that a continuous phase transition accompanies the 
^ I symmetry breaking. The approach used by Eggers was generalized for a larger number of compartments separating 

■ the box. It was shown that in such a case the phase transition becomes discontinuous . 

, Eggers's model is however quite complicated and its analysis must be supplemented with some simplifying assump- 

■ tions. The analysis of fluctuations or dynamical properties seems to be particularly difficult for this model. A possible 
alternative might be to examine a different, presumably simpler model, which would nevertheless grasp the essence 
of the phenomena. 

' O ' The main goal of the present paper is to define and analyze such a model. Our model is a certain generalization of 
5 . the urn model introduced by Ehrenfest In Ehrcnfest model one has balls which are distributed between two urns 
and certain dynamical rules according to which balls can change their location. Ehrenfest model greatly contributed 
to the understanding of some basic notions of statistical mechanics like the equilibrium or approach to the equilibrium. 
Various modifications of this model were examined 

At first sight it seems natural that Ehrenfest model could provide an approximate description of experiments 
; ^ ' with shaken sand. However, none of its generalizations described in the literature exhibits a spontaneous symmetry 
Ci I breaking |ic| ]. The modification which we propose does exhibit a spontaneous symmetry breaking. Although the 
model is rather simple it exhibits a rich behaviour. In addition to having a line of continuous transitions, the model 
has a tricritical point and a line of discontinuous transitions. The discontinuous transitions are screened by very strong 
metastability. Existence of this, unobservable in fact (due to metastability) transition in our dynamical model requires, 
however, some comments. In equilibrium systems it is the free energy which unambiguously determines the phase 
which is stable for a given set of parameters. For our dynamical model we do not know how to calculate the free energy. 
Nevertheless, certain features enabled us to locate the discontinuous transitions in our model. Namely probability 
distribution of the order parameter exhibits characteristics which are typical to equilibrium systems with discontinuous 
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transitions. Additional evidence of such a transition in our model conies form the analysis of characteristic times, 
which we determine solving numerically the first-passage problem. 

Our paper is organized as follows. In section II we introduce the model and examine its steady-state properties. 
In section III, from the solution of the master equation we analyze the steady-state probability distribution and 
fluctuations in the model. Section IV is devoted to the metastable properties and characteristic times obtained from 
the solution of the first-passage problem. Conclusions and possible extensions of our work are discussed in section V. 



II. MODEL AND ITS STEADY STATE 



Our model is defined as follows: N particles are distributed between two urns A and B and the number of particles 
in each urn is denoted as Na and Nb, respectively {Na + Nb = N). Particles in a given urn are subject to thermal 
fluctuations and the temperature T of the urn depends on the number of particles in it. To mimic the effect that a 
larger number of particles leads to a smaller temperature, we use the following rule: 

T{nA^B)=To + A{l-nA,B), (1) 

where ha^b = Na,b/N is a fraction of a total number of particles in a given urn and Tq and A are positive constants. 
With such a choice, temperatures of urns changes linearly with nA,B from the value Tq to Tq + A. For granular 
systems the relation between temperature and the number of particles is very complicated and depends on certain 
parameters like density of particles or the type of driving lO. For the experiment with separation of sand certain 
arguments suggest that T ~ n'^g would be more suitable [Qp^However, with such a choice the location of the critical 
point changes with N . As expected on general grounds, the observed behaviour for finite TV is only a crossover and a 
well-defined phase transition exists only in the limit N ^ oo. Our choice (0), leads to results which are well-defined 
in the thermodynamic limit N oo. However, it would be interesting to examine other choices too. 

Next, we assume that particles within a given urn obey the standard Maxwell-Boltzmann distribution. It means 
that their distribution changes with the height z (above the bottom of an urn) as 

mgNA,B -mgz 

PA,b(z) = — rexp — , 2 

T(nA,Bj T{nA,B) 



where m is the mass of_a particle and g is the Earth acceleration. One can easily calculate that the fraction of particles 

—mgh ] 



obeying distribution (^, which are above a certain height h is given as exp[ ^^^^\ ]. Next, we define dynamics of our 



model as follows: 

(i) One of the N particles is selected randomly. 

(ii) With probability exp[ j^^J''^^ ] the selected particle changes urns, {ua.b is the fraction of particles present in the 
urn of the selected particle). 

In the following we assume that temperature is measured in units such that mgh = 1. It means that our model is 
parametrized only by Tq and A. The Ehrenfest model is recovered in the limit of infinite temperature where every 
selected particle changes urns. 

We define e as a measure of the difference in the occupancy of the urns 

Na~Nb 

For our model, one can find a simple, but exact equation which determines the time average < e > in the steady-state. 
Indeed, in the steady state the flux of particles changing their positions from A to B equals to the flux from B to A. 
It is easy to realize that since the selected particles are uncorrelated, the above requirement can be written as: 

< Na > exp[— -, -1 =< Nb > exp[— —, -1. (4) 

^^T{<Na/N>)^ ^^T{<Nb/N>)^ ^' 

Using the parameter < e > and the relation Na + Nb = N, eq. (^) can be written as 

(i+ < e >)exp[— ; -1 = < e >)exp[— ; -1. (5) 

Of course, eq. admits symmetric solution < e >= for any Tq and A. However, a physically satisfactory solution 
must be stable with respect to fluctuations. To make the stability analysis of e = solution of eq. (p|) wc have to 
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expand it in powers of e. Equating linear terms we obtain that the critical line which separates stable and unstable 
regions for symmetric solution is given as: 

T^-yA72-A/2 (6) 

This line is plotted as a solid line in Fig. [l]. 
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FIG. 1. The phase diagram of the urn model. Regions I and II are symmetric (< e >= 0) and asymmetric (< e >^ 0) 
phases, respectively. In region III (IV) the symmetric (asymmetric) phase is metastable (see text). The solid line is the critical 
line (^. The thick dotted line is the line of first-order phase transitions. The tricritical point (|, ) is denoted as □. 

What is more interesting, however, eq. (^) has also asymmetric solutions with < e >^ 0. Although such solutions 
cannot be written in a closed form, they can be easily determined numerically. An example of such a solution is shown 
in Fig. ^ In addition, expanding eq. (|^) in powers of < e >, one can easily check that in the vicinity of the critical 
point I < e > I (To — 2^o)^^^! which recovers an equilibrium mean- field exponent /? = 1/2 (we assume that A is kept 
const). 

0.6 
0.5 
0.4 

"S 03 

V 

0.2 
0.1 


0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3 
To 

FIG. 2. The absolute value of the order parameter | < e > j as a function of To for A — 0.5. Solid line corresponds to the 
numerical solution of eq. (^). Monte Carlo simulations were made for A*' — 500 (+) and 5000 (*). 

It turns out, however, that in a certain region of phase diagram both symmetric and asymmetric solutions are stable 
(see Fig. ^ It suggests that in that region a line of discontinuous transition should exist. This line would meet a critical 
line at a certain point which, by analogy with equilibrium systems pd| , can be called a tricritical point. Location 
of this point can be easily determined from eq. (p|). Indeed, requiring that the third order (in e) coefRcients are the 

same [|2| and taking in account eq. (||), we obtain the the tricritical point is located at A = |, Tq = . However, 
it is by no means obvious how to determine the line of discontinuous transitions in our model. In equilibrium systems, 
one can calculate the free energy of each phase and this line is determined from the condition that the free energy of 
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both phases are equal. Our model is not an equilibrium system and this procedure cannot be used. Nevertheless, as 
we will show below, behaviour of certain quantities shows striking similarity to equilibrium counterparts. With this 
observation we will be able to locate discontinuous transitions in our model. To address this problem, we have to go 
beyond the steady-state eq. (||), which is the subject of the next section. 

III. MASTER EQUATION 

For our urn model one can write evolution equations for the probability distribution p(Af , t) that in a given urn 
(say A) at the time t there are M particles. These equations easily follow from the dynamical rules: 



piM,t + l) 



p(0,i+l) 



N - M + 1 , , , , Af + 1 , 

-p{M - 1, t)uj{N - M + 1) + — — p(M + 1, t)uj{M + 1) + 



N 



N 



p{M, t){^[l- u:{M)] + ^^^^ [1 - u:{N - M)]} for M = 1, 2 . . . iV - 1 



1 

TV 
1 



p(l,tV(l)+p(0,t)[l-^(iV)], 



p{N, t + l) = -p{N - 1, t)Lu{l) + p{N, t) [1 - Lo{N)], 



(7) 



where w(M) = exp[ypyy^]. For uj{M) = 1 the above equation are equivalent to the ones of the Ehrcnfcst model [||. 

Supplementing the above equations with initial conditions one can solve them numerically [p^ . Examples of such 
solutions in a long-time limit are shown in Fig. |^. 





A=l.35, T|,=0.2 



FIG. 3. Probability distribution p[N{e + obtained as a numerical solution of eq. (j^ in the long time limit (A'^ = 200). 
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As expected, one finds a single- and double-peak distribution in regions I and II, respectively. The width of these 
peaks decreases with the number of particles N. The behaviour of the probability distribution in regions III and 
IV is however more complicated. Although for small N it has three-peak structure, the shape changes in the limit 
N — > oo. The final distribution depends whether the parameters A and To are in region III or IV. In region III, for 
increasing N, the central peak diminishes and in the limit TV — > cxd we obtain a two-peak distribution. On the other 
hand, in region IV only the central peak survives in the thermodynamic limit. As we will show, a line which separates 
these two regions can be interpreted as a line of discontinuous transitions. Before we approximately locate this line, 
let us study fluctuations of the order parameter. Since we know the probability distribution, these fluctuations can 
be calculated numerically, without making additional assumptions Q . Using the variance of the order parameter we 
define the susceptibility k as: 



N < (e- 



1 



N N 
i=0 i=0 



(8) 



where p{..,oo) denotes the steady-state probabilities obtained from eqs. (0). Susceptibility is known to diverge at 
the continuous phase transition. Our calculations for A = 0.5 confirm such a behaviour (Fig. ^) For Tq > 0.25 the 
inverse susceptibility decays linearly which implies 7 = 1 (mean-field value). For Tq < 0.25 (asymmetric phase) 
the behaviour is less clear but we also expect a linear decay of k^^. 

To check our steady-state and master-equation calculations we performed Monte Carlo simulations which are also 
presented in some cases. The implementation of this method for the present model is rather straightforward. 



FIG. 4. The inverse of the susceptibility 




as a function of To for A — 0.5. Solid line corresponds to the numerical 



solution of eq. (^) for A'^ = 5000. Monte Carlo simulations were made for A'^ 
to calculate k we used only a half of the probability distribution. 
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FIG. 5. The scaled susceptibility n/N as a function of A for To — 0.2 obtained from the numerical solution of eq. (|7[) for 
iV = 50, 100, 200, 300, and 500. 

Let us notice that the behaviour of the probability distribution as seen in Fig. ^ bears some similarity to equilibrium 
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systems. In particular, it is known that the simultaneous existence of peaks of both phases is an indicator of the 
first-order phase transitions jlj]. Although our model is defined only using dynamical rules and thus has no free 
energy (at least as defined for equilibrium systems) it does exhibit some other features typical to equilibrium systems. 
To precisely locate the transition point we calculated the scaled susceptibility k/N . We expect that when only a 
single peak survives n/N vanishes, while it remains positive in the double-peak region. For Tq = 0.2 the obtained 
results (Fig. |) locate the transition around A = 1.31. As we will see in the next section there is also a dynamical 
indicator of this transition. 



IV. METASTABILITY AND THE CHARACTERISTIC TIME 



Analysis of the probability distribution for our dynamical model which was done in the previous section shows 
some similarity to equilibrium discontinuous transitions. It is known that such transitions are usually accompanied 
by metastability. In the present section we examine such effects in our model. 

The most transparent indication of metastability is hysteresis. That our model exhibits such a behaviour is clearly 
seen in Fig. |^. One can see, that upon increasing A, the system prepared in the asymmetric state remains in this 
state up to A 1.7 which is close to the limits of stability of this state (see Fig. 0). On the other hand, when 
A decreases the symmetric state persists up to A ~ 1.1 which is again close to the limits of stability of symmetric 
phase. Simulations show that the range of the hysteresis only slightly depends on the system size and the simulation 
time. For large system sizes the model will persist in the initial state until the limits of stability of this state are 
reached, as calculated from the steady-state equation (^ . Such a behaviour indicates a very strong metastability. As 
we will show below, in regions III and IV of the phase diagram certain characteristic times diverge exponentially with 
the system size which corresponds to the broken ergodicity (i.e., dynamical phase space of the model is decomposed 
into disconnected sectors). Let us notice that in (short-ranged) equilibrium systems metastability has always a finite 
lifetime and longer simulations or larger system size shrink the hysteresis. 
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FIG. 6. The absolute value of the order parameter | < e > | as a function of A for To = 0.2 calculated using Monte 
Carlo simulations for N — 2000. Two runs were made with A increasing {+) and decreasing (□), respectively. Due to strong 
metastability no sign of the transition exists for A — 1.31. 



To examine the metastable properties of our model further we will calculate certain characteristic times. It is known 
that for urn models some of these quantities can be calculated as the so-called first-passage problem. We will use a 
similar approach. More specifically, let us consider a configuration where there are M and N — M balls in urn A and 
B, respectively (A^-even). Now, we define an average time t{M,N — AI) needed for such a configuration to reach a 
symmetric configuration (M = N/2). From the dynamical rules of our model one finds that r(M, A^ — Af )'s obey the 
following A. Hnear equations (to simplify notation we omit the second argument in r's): 

t{N) = uj{N)[t{N - 1) + 1] + [1 - uj{N)][t{N) + 1] 
riN - 1) = ^u;{N - l)[r(^ - 2) + 1] + ^L^{l)[riN) + 1] + [1 - ^c.(iV - 1) - ^cviDMN - 1) + 1] 
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-{M)^^c.{M)[t{M-1 



M 
TV 



w(M) 



) + l] 
N - M 



N -M 
N 



uj{N ~ M)[t{M + + 



N 



uj{N ~ M)][t{M + 1)^1] 



t{0.5N + 1) = (0.5 + ^)w(0.5iV + 1) + (0.5 - ^V(0.5A^ - l)[r(0.5^' + 2) + 1] + 

[1 - (0.5 + ^)w(0.57V + 1) - (0.5 - j^)Lj{0.5N - 1)][t(0.57V + 1) + 1] 



(9) 



Similar equations can be written for the Ehrenfest model . 

We were unable to write an explicit solution of the above set of equations. However, its tridiagonal structure 
greatly simplifies its numerical solution. In particular the so-called Gaussian elimination can be used to calculate r's 
recursively [ p^ . Since implementation of this algorithm is straightforward we present only final results. Our unit of 
time corresponds to a single on average update per ball, i.e., we divide the time r calculated from eqs. (^) by N. 

In Fig. ^ we present the A-dependence for t{N). This is the lifetime of totally asymmetric solution. Its rapid 
increase (as a function of TV) in certain range of A indicates the stability of the asymmetric solution. 
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FIG. 7. The characteristic time t{N) as a function of A for (from bottom) = 50, 100, 200, and 300. 
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FIG. 8. The characteristic time t{N) as a function of iV for To = 0.2 and (from left) A = 1.65, 1.7, 1.72, 1.735, 1.74, 
Ae(0.2) = 1.7440675 . . ., 1.77, and 1.8. 



To provide more details on the stability of this solution we have to look at the size dependence of the characteristic 
time and the results are shown in Fig. B-pl. One can see (Fig. ||) that the limits of stability of asymmetric solution. 
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which for Tq = 0.2 equals A = Ac(0.2) = 1.7440675 . . ., separates two regimes. For A > Ac(0.2) t{N) approa ches a 
finite value for increasing N, while it diverges for A < Ac(0.2). From Fig. |9|one can conclude that for A < Ac(0.2), 
t{N) diverges exponentially fast with N. It confirms our previous observations, based on Monte Carlo simulations, 
that a given phase persists up to the limits of its stability. 





14 


/ 7 / 






12 




/ To=0.2 




10 


: 1/ 




z 


8 


1 / / 








1 / / 




o 




-\ / / 




_o 


6 








4 








2 






















2000 4000 6000 8000 10000 12000 

N 

FIG. 9. The characteristic time t(N) as a function of TV for To = 0.2 and (from left) A = 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, and 1.8. 

Form Fig. |o| it follows that at the limits of stability r(A^) increases but as a power of N . From the slope of these 
data we estimate that both for Tq = 0.2 and 0.15 r ^ N^^^ . A different exponent governs the increase of t{N) along 
the critical line (^) In this case we find t{N) ~ 7V^/^, while at the tricritical point N'^/'^ increase is observed (see 
Fig. [l^). It would be desirable to explain such simple power laws using analytical arguments. 
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FIG. 10. The characteristic time t{N) as a function of A'' calculated at (A, To) equal to (solid lines from top): 
(0.25,v'0.25/2 - 0.25/2), (0.5,0.25), (1.7440675,0.2), and (6.775493. .. ,0.15) For two upper lines the points (A.To) are on 
the critical line (^), while for bottom ones on the lines of the limits of stability of the asymmetric solution. The dotted line 
corresponds to the tricritical point. 

The above results concern the behaviour of t{N). We have shown that this quantity can be used to locate the 
critical line or the limits of stability of the asymmetric solution. However, this quantity docs not provide any indication 
of the discontinuous phase transition which we located between regions III and IV on Fig. |l|. To have a dynamical 
indication of this transition we have to look for other characteristic times. First, let us notice (Fig. pd| ) that in region 
II the distribution of t's as function of M is rather flat which means that there is basically a single time scale in the 
model. In regions III and IV, however, one can see a more pronounced variability of t{M). In Fig. |l^ we show that 
t(0.57V + 1) (which is the shortest of our t's) can be used to locate a discontinuous phase transition. Indeed, this 
quantity increases exponentially fast with N for A < 1.31 and remains bounded for larger values of A. The change 
of behaviour of t(0.5A^ + 1) coincides, within our numerical accuracy, with the change of the probability distribution 
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as shown in Fig. || 
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FIG. 11. The distribution of characteristic times t{M) as a function of M calculated To — 0.2 and (from top) A = 0.9, 1.2, 
1.5. 
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FIG. 12. The characteristic time r(0.5A'' + 1) as a function of N calculated To = 0.2 and (from top) A = 1.27, 1.29, 1.30, 
1.31, and 1.33. 

Let us also notice that the time scale reached in the numerical solution of eqs. (||) ('^ 10^"') is much larger than 
what is accessible using Monte Carlo simulations. The main limitation in our calculations is finite accuracy of the 
computations. Using longer representations of real numbers (we used the FORTRAN real*8 type ) one can study 
much larger time scales. 

V. CONCLUSIONS 

In the present paper we examined an urn model which undergoes a symmetry breaking transition. Our work was 
motivated by recent experiments on the separation of shaken sand. The proposed model exhibits a rich behaviour 
(continuous/discontinuous transitions, a tricritical point, metastability) and its properties could be reliably deter- 
mined. There are experimental indications that in the three-compartment case the separation of shaken sand is a 
discontinuous transition [Q. However, due to strong metastability only the hysteretic behaviour is observed, similarly 
as in our simulations in Fig. |6[ Our work suggests, that there are some indicators (probability distributions, charac- 
teristic times) which could be used to locate a discontinuous transition in such systems. It would be interesting to 
look for such indicators in experimental systems. 

As an extension of our work, one can examine models with a different relation between the effective temperature and 
the number of particles or models with larger number of compartments. Even further extensions are related with a 





9 



much different interpretation of our model. One can imagine that two compartments A and B are two groups of people 
(e.g., supporters of certain presidential candidates). It seems reasonable to assume that under certain conditions a 
larger number of people in a given group will increase attractive bonds within this group (energy). On the other hand 
it will also increase a probability that a certain member will leave the group (entropy). Competition of this two effects 
is likely to produce a similar symmetry breaking transition to the one described in the present paper. However, since 
people are known to create various local connections (e.g. within a group) it would be desirable to take such effects 
into account too. In our model particles are chosen randomly and such correlations are omitted. 

ACKNOWLEDGMENTS 

This work was partially supported by the Swiss National Science Foundation and the project OFES 00-0578 
"COSYC OF SENS". 



[1] P. B. Umbanhowar, F. Melo, and H. L. Swinney, Nature 382, 793 (1996). 
[2] T. Shinbrot and F. J. Muzzio, Nature 410, 251 (2001). . 

[3] H. J. Schlichting and V. Nordmeier, Math. Naturwissenschaften Unterr. 49, 323 (1996) (in German). 

[4] J. Egger, Phys. Rev. Lett. 83, 5322 (1999). 

[5] V. Kumaran, Phys. Rev. E 57, 5660 (1998). 

[6] I. Goldhirsch and G. Zanetti, Phys. Rev. Lett. 70, 1619 (1993). 

[7] D. van der Meer, K. van der Weele, and D. Lohse, Phys. Rev. E 63 061304 (2001). K. van der Weele, D. van der Meer, 

M. Versluis, and D. Lohse, Europhys. Lett. 53, 328 (2001). 
[8] P. Ehrenfest and T. Ehrenfest The Conceptual Foundations of the Statistical Approach in Mechanics (Dover, New York, 

1990) M. Kac and J. Logan, in Fluctuation Phenomena ed. E. W. Montroll and J. L. Lebowitz (North-Holland, Amsterdam, 

1987). 

[9] see e.g., R. Metzler, W. Kinzel, and I. Kanter, J. Phys. A 34, 317 (2001); C. Godreche and J. M. Luck, J. Phys. A 32 
6033 (1999) and references therein. 
[10] Certain urn models do exhibit a symmetry-breaking transition. However, their essential ingredient is an infinite number of 
urns which is not suitable for our purposes. For details see: P. Bialas, L. Bogacz, Z. Burda and D. Johnston, Nucl. Phys. B 
575, 599 (2000). 

[11] M. Blume, V. J. Emery and R. B. Grifiiths, Phys. Rev. A 4, 1071 (1971). 

[12] One can easily notice that even-order coefficients in the expansion of both sides of eq. ^ are the same and thus cancel 
each other. 

[13] Analytical approach to evolution equations (^) seems to be very difficult since the variable M appears in an exponential 
factor ul!{M). 

[14] K. Binder and D. P. Landau, Phys. Rev. B 30, 1477 (1984). 

[15] A. Lipowski, J. Phys. A 30, L91 (1997). K. P. N. Murthy and K. W. Kehr J. Phys. A30, 6671 (1997). 
[16] R. Sedgewick, Algorithms, (Addison- Wesley Publishing Company, New York, 1988). 

[17] Actually, only a right-hand (respective to the tricritical point) portion of the curve (^) can be called a critical line and 
only on this part the scaling t{N) ~ N^^'^ is observed. The left-hand side is only the limit of the stability of the symmetric 
solution. Since the asymmetric solution is stable on this part, t{N) diverges exponentially with N . 



10 



